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Dissolved ions can alter the local permittivity of water, nevertheless most theories and simulations 
ignore this fact. We present a novel algorithm for treating spatial and temporal variations in the 
permittivity and use it to measure the equivalent conductivity of a salt-free polyelectrolyte solution. 

Our new approach quantitatively reproduces experimental results unlike simulations with a constant 
permittivity that even qualitatively fail to describe the data. We can relate this success to a change 
in the ion distribution close to the polymer due to the built-up of a permittivity gradient. 


The dielectric permittivity 5 measures the polarizabil¬ 
ity of a medium subjected to an electric field and is one of 
only two fundamental constants in Maxwell’s equations. 
The relative permittivity of pure water at room temper¬ 
ature is roughly 78.5, but charged objects dissolved in 
the fluid significantly reduce the local dielectric constant 
because water dipoles align with the local electric field 
created by the object rather than the external field 
When ions accumulate in the vicinity of a charged object 
they further reduce the local dielectric constant Emia, 
causing a permittivity gradient that repels ions from the 
surface ISHS]. The screening of electrostatic forces be¬ 
tween charged objects is therefore affected, and it also 
influences any properties that depend on the specific ion 
distribution. However, almost all computational and the¬ 
oretical work to date using an implicit water model as¬ 
sumes a constant dielectric permittivity. This is partially 
due to a lack of suitable numerical approaches, since only 
few electrostatic algorithms can include spatial changes 
in the dielectric permittivity [inHi3]. 

Electrophoresis is the directed motion of an object in 
an aqueous solution subject to an external electric field. 
The relative ease with which electric fields can be applied 
experimentally has led to wide use of electrophoresis in 
the characterization of polymers pTUIB] , colloids [nHii], 
and cells [22j [23l , which tend to ionize in aqueous solu¬ 
tions. Measuring the electrophoretic velocity of individ¬ 
ual particles is often difficult from a technical standpoint. 
For this reason, polyelectrolyte solutions are often char¬ 
acterized by their conductivity [241433] . 

In both electrophoresis and conductivity, the distribu¬ 
tion of oppositely charged counterions around the object 
determines the magnitude of the relative velocity between 
the object and the fluid. The surrounding counterion 
cloud is comprised of two separate layers: the Stern layer 
and the Debye layer [34]. The Stern layer, often called 
the stagnant layer, consists of strongly adsorbed ions ad¬ 
jacent to the charged object that reduce the effective sur¬ 
face charge. Beyond the Stern layer is the Debye layer, 
also called the diffuse layer, which consists of ions that 
are free to move relative to the charged surface. The two 
layers are collectively referred to as the electric double 
layer (EDL). 

For the study of charged macromolecules a coarse¬ 


grained approach is necessary due to the excessive sys¬ 
tem size. Molecular Dynamics (MD) simulations of 
charged systems typically work with the restricted prim¬ 
itive model by simulating the ions as hard spheres while 
accounting for the solvent implicitly through a con¬ 
stant background dielectric. Crucially, the solvent me¬ 
diates hydrodynamic interactions and reduces the elec¬ 
trostatic interactions due to its polarizability. Hydro- 
dynamic interactions significantly impact the electroki- 
netic properties and lead to qualitatively different be¬ 
havior [35l [36]. Polarizability actually depends on the 
local electric field [371138]: these local variations in the 
polarizability thicken the Debye layer by repelling the 
counterions (MSI HSlEnillOI. More detailed atomistic 
computer simulations and experiments have also indi¬ 
cated pronounced deviations to standard theories for ion 
distributions around charged objects if the continuum 
solvent approach is replaced by an explicit water envi¬ 
ronment mnHini. Hence, these findings show that the 
molecular properties of water are important for a detailed 
description of the EDL. However, the consequences for 
electrokinetic properties are largely unknown. 

In this letter, we show that including spatially and 
temporally varying dielectric properties significantly in¬ 
fluences the structure of the Debye layer and the dynamic 
properties of polyelectrolytes. By means of a novel algo¬ 
rithm, we locally couple the permittivity to the local ion 
concentration. When applied to a polyelectrolyte solu¬ 
tion in an external electric field, the algorithm quantita¬ 
tively reproduces experimental data for the conductivity, 
while simulations assuming a constant dielectric back¬ 
ground disagree qualitatively with experiment. We show 
that the decreased permittivity in the vicinity of the poly¬ 
electrolyte reduces the fraction of condensed counterions 
on the polyelectrolyte backbone, explaining the experi¬ 
mentally observed increase in the equivalent conductivity 
at high monomer concentrations. 

We performed standard coarse-grained simulations us¬ 
ing a Weeks-Chandler-Anderson (WCA) potential [50| 
for steric interactions with an equilibrium distance of 
a = 0.3 nm between particles. Adjacent monomers 
are connected by nitely extensible nonlinear elastic 
(FENE) bonds. Simulations looking at static ion dis¬ 
tributions were performed with a Langevin thermostat. 
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Figure 1. (color online) Our system setup with polyelec¬ 
trolytes (black spheres) and counterions (red spheres). The 
scheme to calculate local charge concentrations is depicted 
by the lattice in the top left corner. The concentration of 
a lattice cell is determined via a weighted summation of all 
charged particles in the 7^ surrounding MEMD lattice cells. 
Ions thus influence the local permittivity within a distance 
d = 1.4 nm, or two Bjerrum lengths. 


while simulations in which dynamic quantities were mea¬ 
sured used the D3Q19 lattice-Boltzmann (LB) fluid in 
ESPResSo [5TH53] and the corresponding thermostat. A 
detailed description of the simulation method can be 
found in the supplemental material [54]. 

The interactions between charges are calculated us¬ 
ing an extension of a local electrostatics algorithm that 
was first introduced by A. Maggs [55] and later inde¬ 
pendently adopted for MD simulations by J. Rottler [56| 
and 1. Pasichnyk m- This method of calculating elec¬ 
trostatic interactions has been coined Maxwell Equations 
Molecular Dynamics (MEMD). This algorithm’s locality 
permits arbitrary changes in the dielectric constant [58| . 
The motion of charges qiVi is interpolated to a regu¬ 
lar lattice. The resulting electric current gives rise to 
a change in the magnetic field B and the displacement 
field D = eE, following Maxwell’s equations. The algo¬ 
rithm, however, treats the propagation speed of the mag¬ 
netic field c as a tunable parameter. For a wide range of 
values of c, MEMD produces the correct particle dynam¬ 
ics and statistic observables ISi Eg. The algorithm is 
therefore computationally efficient, since only two update 
steps for the electromagnetic fields are performed per MD 
time step. The permittivity 5 is used twice: The coupling 
of the displacement field to the magnetic field and vice 
versa (Ampere’s law and Faraday’s law). Effectively, this 
creates two new driving forces if 5 depends on space and 
time: The influence of magnetic waves that are partially 
reflected in lattice cells with variable permittivity, and 
a force directly pointing in direction of the permittivity 
gradient, following the Lorentz force Fl = qD/e-\-v x B, 

The MEMD grid spacing is aMEMD = 0.4 nm and the 
time step is set to A^memd = A^md- The last param¬ 
eter is an artificial mass /mass = 0.05mo. To take into 
account changes in the solvent polarization due to ions 


in the double layer laiiE], the local dielectric permit¬ 
tivity is calculated from the nearby ion density using the 
empirical function found by B. Hess et al. |4], 


1 + 0.278-C’ ^ ^ 

where C is the molar salt concentration in [mol/L] or 
[M]. For each lattice cell, the charge density is averaged 
in the surrounding 7^ cells as shown in the top left of 
Figure weighted by the inverse square of the depicted 
shell number, resulting in weights from 1 to 1/16. This 
is because the electric field created by an ion decays as 
1/r^ with the distance r. Assuming linear response, the 
polarizability should be proportional to this local electric 
field. The weighting also guarantees that charges enter¬ 
ing or leaving the cube do not lead to sharp jumps in the 
dielectric constant, which cause unphysical behavior such 
as jumps in the total electric field energy. The volume 
taken into account for the averaging is (2.8 nm)^, which 
for a Bjerrum length of /p ~ 0.7 nm is roughly the extent 
over which the electrostatic interactions are significant 
compared to thermal fluctuations. 

The temporal changes in permittivity lead to addi¬ 
tional interactions which have been called spurious m, 
but are indeed physical as pointed out by Rottler and 
Maggs m and Pasichnyk et al. m- In our simulations, 
however, these effects are very small since the algorithm 
described above only allows smooth and slow changes in 
the local permittivity. 

Initially we examine the counterion distribution 
around a single rod-like chain of = 80 monomers sepa¬ 
rated by 0.3 nm in a cubic box with a side length of 24 nm 
with periodic boundary conditions. In our conductivity 
simulations, we vary the monomer concentrations C by 
placing M G {1,..., 25} chains of length N = 30, 45, and 
60 in a box of size (32 nm)^, depending on the desired 
polyelectrolyte concentration. On a modern workstation 
(Intel i7-5820K CPU, Nvidia GTX 780 Ti GPU), the 
simulation time for each curve in Figure 4 was around 46 
hours. This includes 16 points with 3 different polymer 
lengths for each, adding up to a total of 48 simulation 
runs during this time. 

We first verify the validity of our new approach, in 
which the dielectric constant dynamically adapts to the 
local ion concentration. To this end, we use the itera¬ 
tive scheme sketched in Figure to calculate the coun¬ 
terion distribution around an infinite rod-like polyelec¬ 
trolyte fixed in space. 

We start out with a constant permittivity of 5 = 78.5 
and obtain a counterion distribution from an equilibrium 
MD simulation. We then map the cylindrically symmet¬ 
ric salt concentration to a local permittivity using equa¬ 
tion 0 and used it as a fixed permittivity for the sub¬ 
sequent simulation run. A successive under-relaxation 
scheme equally weighting the two preceding results con¬ 
verges to the counterion distribution in Figure We 
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Figure 2. Iterative scheme: Starting with a constant back¬ 
ground permittivity, the resulting ion concentrations of suc¬ 
cessive simulations are mapped to a spatially varying but hxed 
permittivity distribution, until the scheme converges towards 
a stable equilibrium. 


then apply the new time-dependent adaptive scheme, 
seen in Figure to calculate the local charge concentra¬ 
tion and permittivity to a fixed rod-like polyelectrolyte, 
as well as a fully flexible polyelectrolyte. Figure dis¬ 
plays all three counterion distributions. 
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Figure 3. (color online) Counterion distribution around the 
polyelectrolyte backbone with varying permittivity as a func¬ 
tion of the distance from the polyelectrolyte surface. The iter¬ 
ative scheme for a stiff rod (dashed line), the adaptive scheme 
for a stiff rod (dotted line), and the adaptive scheme for a flex¬ 
ible polymer (red triangles) are almost identical. They differ 
qualitatively from the analytical Poisson-Boltzmann (PB) so¬ 
lution (thick black line) for a uniform permittivity £ = 78.5. 
The corresponding permittivity £ (inset) in the iterative case 
goes from 78.5 in the bulk to 60 close to the polyelectrolyte 
backbone, similar to what is observed near a charged sur¬ 
face [Wi Pi- The counterion distribution for a sharp rise in 
the permittivity from 2 within the rod to 78.5 in the fluid 
(thin light green line) shows only minor deviation from the 
PB solution. 

The non-monotonicity in the counterion distribution 
is not predicted by Poisson-Boltzmann (PB) theory with 
a fixed constant dielectric constant. Note that in the 
case of a uniform background permittivity there is also 
a sharp increase in the counterion density near the poly¬ 
electrolyte backbone, however, the increase only ex¬ 
tends over a relatively short distance (approximately 
O.lcr = 0.03 nm, data not shown) compared to the simu¬ 
lations with a varying dielectric constant (approximately 
4(7 = 1.2 nm). 


The physical basis of the extended initial increase in 
the counterion concentration when the permittivity is 
adapted locally is the permittivity gradient in the vicin¬ 
ity of the polyelectrolyte, where there is an increase from 
approximately 5 = 60 at the polymer surface to 5 = 78.5 
in the bulk (see inset of Figure]^. The observed distribu¬ 
tion is the combined result of the repulsion of counterions 
by the permittivity gradient, the electrostatic attraction 
of the ions to the polyelectrolyte backbone, and ther¬ 
mal fluctuations. Interestingly, there is almost no visi¬ 
ble difference in the counterion distribution around the 
flexible polyelectrolyte and the infinite stiff charged rod, 
demonstrating that the cell model is a good approxima¬ 
tion [63l l64) . 

Our results resemble earlier observations for the coun¬ 
terion distribution around a colloid with spatially varying 
dielectric background laizi- All simulations with varying 
permittivity also agree qualitatively with atomistic sim¬ 
ulations of a similar system using a Kirkwood-Buff based 
force field for NaCl [65] in combination with the extended 
simple point charge model (SPC/E) water model [66] . 
The atomistic simulations displayed a similar depletion 
of counterions very close to the polyelectrolyte backbone 
and a subsequent rise of the distribution. The excellent 
agreement between our three simulation setups, and both 
existing ion distributions for colloids and our atomistic 
simulations, demonstrates that our method for dynami¬ 
cally adapting the local permittivity produces physically 
very reasonable results. 

In Figure we have also plotted simulation results 
where there is a sharp dielectric interface (thin green line) 
at the surface of the polyelectrolyte, with a discrete rise 
from 2 within the polyelectrolyte to 78.5 in the fluid. 
As observed in other studies [Sin], the discrete change 
in the dielectric interface also produces a thickening of 
the Debye layer. However, the difference in the distri¬ 
bution compared to the Poisson-Boltzmann (thick black 
line) result is significantly less than in our simulations 
with a permittivity adapted to the local salt concentra¬ 
tion (dotted line). This shows that one not only needs 
take into account the reduced permittivity within a poly¬ 
mer, colloid, or charged surface, but also the reduction 
in the permittivity in the surrounding Debye layer, i. e. a 
gradient in the permittivity. 

The structural differences within the EDL in Eigurej^ 
have little influence on many properties such as the ra¬ 
dius of gyration or the polymer diffusion coefficient m- 
However, we found that adapting the local dielectric con¬ 
stant significantly impacts the response of the system to 
an external electric field. This is because the electroki- 
netic behavior of the system strongly depends on the hy¬ 
drodynamic and electrostatic friction between the poly¬ 
electrolyte backbone and its counterions, and is thus very 
susceptible to changes within the EDL. 

We simulated the equivalent conductivity A (the con¬ 
ductivity over the polyelectrolyte concentration) of a 
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Figure 4. (color online) The rescaled equivalent conductiv¬ 
ity A/Ao over monomer concentration \/C for simulations 
with constant permittivity (green dashed, blue dotted) and 
locally varying (red solid) permittivity. Experimental data 
(gray symbols) from Kwak and Hayes [28], Colby et ah [^ . 
and Lipar-Ostir et ah m is reproduced with locally vary¬ 
ing permittivity, while we observe a qualitative difference for 
a constant dielectric background. The reason for this is a 
significant drop in the dielectric constant around the poly¬ 
electrolyte backbone, as seen in Figure [^ and the subsequent 
repulsion of counterions from the polymer. 

salt-free solution of polyelectrolytes as a function of the 
monomer concentration C as plotted in Figure [^together 
with experimental data. The equivalent conductivity is 
simply: 


where j is the current density, E is the applied electric 
field, and C is the monomer concentration. The con¬ 
ductivities have been additionally normalized by an ex¬ 
trapolated Ao((7 = 0), since the hydrodynamic radius of 
the specific ions strongly influences the diffusion and con¬ 
ductivity even at infinite dilution. We performed two sets 
of simulations with a constant background permittivity, 
one with Sr = 78.5 and one where the ion concentration 
in the simulation box was substituted into equation 
to determine the background permittivity. Both sets of 
simulations with a constant dielectric background show 
a monotonic decrease in the equivalent conductivity, in 
good agreement with scaling theories which ignore the 
effect of dielectric contrast [24l [25l [271 EH [3Ql [68] . More 
importantly, our simulations that adapt the permittivity 
to the local ion concentration quantitatively reproduce 
the unexpected rise in conductivity at high salt concen¬ 
trations. 

From the snapshots in Figure at monomer concen¬ 
trations of 1 mM, 10 mM, and 100 mM, we see a coiling of 
the polyelectrolytes (black spheres) as well as a decrease 


in the fraction of condensed counterions (red spheres) 
at high concentrations for the simulations including di¬ 
electric variations (top right snapshot). The simulations 
assuming a constant dielectric background (bottom right 
snapshot) resulted in a larger number of condensed coun¬ 
terions (blue spheres). To quantify the fraction of con¬ 
densed counterions, we used the criterion suggested by 
L. Belloni [69] and M. Deserno IZO]. In Figure]^ we plot 
the fraction of condensed counterions fed as a function 
of the monomer concentration C. 
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Figure 5. (color online) The local permittivity £poiy around 
the polymer backbone (black squares), and the fraction of 
condensed counterions fed for constant dielectric background 
(blue circles) and varying local permittivity (red triangles) 
as a function of the monomer concentration. The maximum 
in the fraction of condensed counterions is due to the sharp 
decrease in the permittivity with increasing monomer concen¬ 
tration and occurs at the same concentration as the minimum 
in the equivalent conductivity in Figure [^ 

While fed monotonically increases for simulations with 
a constant dielectric background, we observe a maximum 
and subsequent decrease when we adapt the dielectric 
constant to the local ion concentration. The maximum 
occurs at the same monomer concentration as the min¬ 
imum in conductivity in Figure The decrease in feei 
results in a higher effective charge of the polyelectrolytes, 
which both increases their mobility and the net charge 
they carry with them. In addition, there is an increase in 
the number of free counterions contributing to the overall 
conductivity. This is why the maximum in the fraction 
of condensed counterions results in a minimum in the 
equivalent conductivity of the solution. 

We average the dielectric permittivity within 1.4 nm 
of the polyelectrolyte backbone and plot the results as 
black squares in Figure We find that the reason for 
the drop in feei is a decrease in dielectric permittivity and 
therefore a steeper gradient in 5, resulting in a stronger 
dielectric repulsion from the backbone. The decrease in 
5 is enhanced due to the coiling of the chains at high 
polyelectrolyte concentrations as seen in the snapshots 
in Figure This increases the density of ions within 
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the polymer coil, which leads to the large decrease in the 
local permittivity around the polyelectrolyte backbone. 

The dependence of the conductivity on the monomer 
concentration is largely independent of the specific sys¬ 
tem being studied. This can be seen by the compari¬ 
son of three different experimental data sets in Figure]^ 
Colby et al. m measured the equivalent conductiv¬ 
ity poly(styrenesulfonate) (PSS) with Na+ counterions 
at 296 K with a molecular weight of M = 1.2 x 10^. In 
contrast, Kwak and Hayes [28] used poly(styrenesulfonic 
acid) (PSA) with Li+, Na+, K+, Cs+ at 298 K with 
a molecular weight of M = 5.6 x 10^. Finally, Lipar- 
Ostir et al. [26] measured the equivalent conductivity 
of poly(anetholesulfonic acid) (PAS) with H+ measuring 
seven sets at temperatures between 278 and 308 K. De¬ 
spite the different counterion types, temperatures, poly¬ 
electrolytes, and polyelectrolyte lengths, all experimental 
data sets collapse onto our simulations that include vari¬ 
ations of the local dielectric constant. This demonstrates 
that our model captures the essential physics. 

Our results demonstrate that the inclusion of local 
variations in the dielectric permittivity in coarse-grained 
MD simulations strongly influences the static and dy¬ 
namic properties of charged macromolecules in aqueous 
solutions and should be taken into account in the mod¬ 
eling process. We expect that dielectric contrast is a key 
factor in determining the electrokinetic properties of soft 
matter systems, such as colloid and polyelectrolyte elec¬ 
trophoresis; electroosmotic flow; and induced charge elec¬ 
troosmosis and electrophoresis, and that our algorithm 
can be applied successfully to all these systems. 
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